Ratio control in a cascade model of cell differentiation 
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We propose a kind of reaction-diffusion equations for cell differentiation, which exhibits the Turing 
instability. If the diffusivity of some variables is set to be infinity, we get coupled competitive 
reaction-diffusion equations with a global feedback term. The size ratio of each cell type is controlled 
by a system parameter in the model. Finally, we extend the model to a cascade model of cell 
differentiation. A hierarchical spatial structure appears as a result of the cell differentiation. The 
, size ratio of each cell type is also controlled by the system parameter. 
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Reaction-diffusion systems are important models for various nonlinear processes in nonequilibrium physics and 
biological systems Turing first proposed a reaction-diffusion model to explain pattern formation in developmental 
biology ■ Various pattern formation has been studied in theoretical models [1, [i| . 

There is a characteristic wavelength in the Turing pattern. Kondo and Asai observed that the number of stripes in 
the skin pattern of a tropical fish increases to keep the characteristic wavelength constant, when the fish grows with 
time In the skin pattern, the periodic pattern is constructed of pigment cells with different colors. The different 
types of pigment cells are mutually competitive. 

On the other hand, there are several observations that various types of cells are differentiated during the develop- 
■ mental process, but the number ratio of the different cell types does not change very much as the body size increases. 
q-i A typical example is the cellular slime mold Dictyostelium discoideum. The amoebic state changes to the fruiting 
'— body via the migrating slug state, if the breeding condition becomes worse. Prespore cells and prestalk cells are 
differentiated in the process. The number ratio is kept almost constant even if the body size is changed [6|]. The 
ratio control and the pattern formation are performed in two steps. In an early stage, the ratio of the two-type cells 
j^. ■ is regulated by DIF (Differentiation Inducing Factor), and later prestalk cells move into a tip region within the slug 
{Nl | through a chemotaxis by the cAMP 0, HI]. The two cell types, i.e., the prespore and the prestalk are competitive 
m this system. Another example is the differentiation of blood cells from the stem cells in the bone marrow. The 
number ratio of the red blood cells, the white blood cells, and the platelets is kept to be roughly constant. 
t^J- i Another typical rule of the cell differentiation is a cascade control of the differentiation process by complicated gene 
networks. In the segmentation process of the Drosophila, the differentiation proceeds from a large scale to a small 
ON i scale in a hierarchical manner. A cascade network of genes and proteins such as bicoid protein, gap genes and pair-rule 
genes are identified in detail @. Several competitive relations between two genes are known also in this process. For 
example, the gene engrailed (en) for the posterior compartmental specification and the gene wingless (wg) for the 
anterior compartment are mutually competitive. Several theoretical models of the hierarchical gene network were 
proposed for the segmentation process [l(| EH • 

The analyses of the specific gene network and the pattern formation for each system are important in the devel- 
opmental biology. However, in this paper, we consider a very simple cascade model of competitive reaction-diffusion 
equations and propose a mechanism of the ratio control from a view point of the nonlinear dynamics. 

II. RATIO-CONTROL IN A COMPETITIVE REACTION-DIFFUSION SYSTEM 

If two cell-types are mutually competitive, only one type of cells are locally activated by the lateral inhibition. 
The two cell-types are bistable, and a homogeneous state of one-type cell will appear. However, if some long-range 
enhancement of the competitive type is included in the system, two types of cells will appear in a large system. 
Meinhardt and Gierer proposed a model of the mutual induction of locally exclusive states, and generated a stripe- 
like pattern [l2| . We study first a simplified version of their model. The model is written as coupled reaction-diffusion 
equations: 

dX 1 cXi + dY 2 2 

~df l + aX 1+ bX 2 - Xl+DxV Xl ' 
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FIG. 1: (a) Eigenvalue \k for Eq. (2) at a = 0.4, c = 2,d = 1, flx = 1, Dy = 20 and 6 = 5 (solid curve), 6 = 4 (dashed curve), 
and 6 = 3 (dotted curve), (b) Stationary profiles of X\ (solid curve) and X2 (dashed curve) at a = 0.4, 6 = 5, c = 2,d = 1 and 
D x = l. 
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FIG. 2: (a) Eigenvalue for Eq. (4) at a = 0.4, c = 2, d = 1,-Dx = 1 and b — 5 (solid curve), 6 = 3 (dashed curve), and 
6=1 (dotted curve), (b) Stationary profiles of X\ (solid curve) and X2 (dashed curve) at a = 0.4, 6 = 5, c = 2, d = 1, Dx = 1. 
(c) Time evolution of X\ (solid curve) and X2 (dashed curve) from X\ = 0.4 + 0.04cos(27ra;/Z/) + 0.1 cos(47ra;/L) and X2 = 
0.4 - 0.1cos(4ttx/L). 



^ = CX2 + m -X 2 + D X V>X 2 , 
dt l + aX 2 + bX 1 



dY 1 

~~dt 



X 1 -Y 1 +DyV J Y u 



dY 

= X 2 -Y 2 + D Y V 2 Y 2 , (1) 

where a, 6, c, d are positive parameters, X% and X 2 obey the Hill type equation, Y\ and 1^ are produced respecti vely 
from Xi and X 2 . The Hill type equation is often used in the biochemistry and the gene regulatory network [l^, Il4j . 
The enhancement and the suppression of the reaction are expressed in the numerator and the denominator in the first 
term of on the right-hand side of the Hill equation. That is, both Xi and X 2 acts as activators for themselves by the 
terms of cX\ or cX 2 , and acts as inhibitors for the other by the term of bX 2 and bX\. The two components X\ and 
X 2 compete with each other when b is sufficiently large. The terms aX\ and 0X2 in the denominators express the 
effect of the saturation. We assume that X\ and X 2 correspond to proteins (morphogens) determining competitive 
cell types such as the pigment cells of different colors in the skin pattern, or the prestalk and the prespore in the slime 
mould. The term of dY 2 and dY\ implies that Y 2 and Y\ are activators respectively for X\ and X 2 . We assume that 
Dy ^> Dx- It means that the morphogens Y\ and Y 2 diffuse rapidly and work as long-range activators respectively 
for X 2 and X\. 

There is a uniform solution Xi = X 2 —Y\ — Y 2 = X = (c + d — l)/(a + b) to Eq. (1). It is expected that X 2 is 
depressed when X\ increases locally and vice versa, because X\ and X 2 are mutually competitive. If the perturbation 
of the form: X\ = X a + 8x\ cos(fcx), X 2 = Xq — 6x1 cos(kx), Yi = X + 8y\ cos(kx), Y 2 = X — Syi cos(kx) are assumed 
owing to the competitive relation, the deviations 5x\ and Syi satisfy 

dSxx _ cSxi - dSyi (c + d)(b - a)X Q Sxi 2 
~dF l + (a + b)X + {l + ( a + b)X o y U + Uxk > dXl > 
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FIG. 3: (a) Stationary profiles of X\ (solid curve) and X2 (dashed curve) at a = 0.4, b — 5, c = 2, di = 0.4, d% — 1 and Dx = 1 
for Eq. (6). (b) Size ratio of the Xi-dominant region as a function of di for a — 0.4, b = 5, c = 2, di = 0.4, di — 1 and Dx = 1. 
The rhombi denote numerical results and the dashed curve is di/(di + ^2). 



6 Xl - {1 + D Y k 2 )5 yi . (2) 



The eigenvalue A^ can be calculated from Eq. (2). Figure 1(a) displays A& for 6 = 5 (solid curve), 6 = 4 (dashed curve) 
and 6 = 3 (dotted curve) at a = 0.4, c = 2, d = 1, Z?x = 1 and Dy = 20. (The parameter values a, b, c, d, Z?x and D y 
are not realistic biological ones, but we study the nonlinear system as just a model system in this paper.) The Turing 
type instability is expected for a finite wavenumber. We have performed a one-dimensional numerical simulation. 
The system size L = 100, and the periodic boundary conditions are imposed. Figure 1(b) shows stationary profiles 
of X\{x) (solid curve) and X 2 (x) (dashed curve) at a = 0.4,6 = 5,c = 2,d = l,Dx = 1 and Dy = 20. A spatially 
periodic pattern appears, in which the AVdominant region and ^-dominant region reciprocate spatially. This is a 
Turing pattern in the competitive reaction-diffusion systems. Note that the Y component with large diffusivity acts as 
activators for the competitor, and it facilitates the competitor in spatially distant regions. It is different from the usual 
Turing pattern in usual activator-inhibitor systems, where the inhibitor with large diffusivity inhibits the activator. 
The mechansim of the pattern formation by the long-range activation of the competitor is essentially the same as the 
model proposed by Meinhardt and Gierer almost 30 years ago [r2j . Kondo and Asai performed a numerical simulation 
of the usual type of activator- inhibitor system to explain the skin pattern of the tropical fish [f| , but our model might 
be more plausible for the pattern formation, because the pigment cells are mutually competitive. 

Next, we generalize the model equation (1) to a model equation in which the whole system is separated into two 
regions of Xi-dominant region and ^-dominant region. If Dy is infinitely large and a steady state is obtained for 
the Y component, Eq. (1) is reduced to 



dXi cX 1 + d{X 2 



dt 1 + aXi + bX-, 



Xi +D X V 2 X 1 , 
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dX 2 cX 2 + d(X 1 ) , 

where {X) is the spatial average of X, i.e., (X) = (1/L) Xdx in a one-dimensional system and {X) = 

(1/L 2 ) Jq Xdxdy in a two-dimensional system. It is because Y becomes uniform owing to the infinitely large 
diffusivity of Y in Eq. (1). This is our original model. For this model equation, the small deviation 5x\ from the 
uniform solution satisfies 

dSx! c5x! (c + d)(b-a)X Q Sx 1 2 

~dT = l + (a + 6)A + {l + (a + 6)X P " (1 + Dxk (4) 
for k 7^ 0, because Syi is assumed to be zero owing to the uniformity of Y. On the other hand, for k = 0, 

dSx\ [c — d)8xi (c + d)(b — a)Xo5xi 



dt l + (a + b)X {l + {a + b)X y 



5xi. (5) 



Fugure 2(a) displays the eigenvalues Afc for 6 = 5 (solid curve), 6 = 3 (dashed curve) and 6=1 (dotted curve) at 
a = 0.4, c = 2, d = l,Dx = 1. The eigenvalue Xk at k = is negative for these parameter values, therefore, the 
Fourier mode with k = is stable. The uniform solution is unstable for 6 = 5 and 6 = 3. The eigenvalue takes the 
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largest value for the smallest wave number k = 2tt/L. Figure 2(b) shows a stationary solution for 6 = 5 for the one- 
dimensional system with periodic boundary conditions of size L = 100. The initial condition was X\ = 0.5 — 0.1(x/ L) 
and X 2 — 0.4 + 0.1(x/L). For this initial condition, the whole space is separated into the Xi-dominant region at 
x < L/2 and the ^-dominant region at x > L/2. There are domain walls between the two domains. However, even 
if the initial condition is random, the whole space is separated into the two regions with the same size, although the 
position of the AVdominant region becomes random. It is because Xk takes a maximum value at k = 2it/L. If the 
initial condition is X\ — 0.4 + 0. 04 cos(27nr/L) + 0.1 cos(47ra;/L) and X2 = 0.4 — 0.1 cos(47ra;/L), four domains, i.e., two 
Xi-dominant regions and two ^-dominant regions appear initially, however, there is an attractive interaction between 
the domain walls and finally the two-domain structure is obtained as shown in Fig. 2(c). The time evolution is similar 
to the coarsening process in the one-dimensional Ginzburg-Landau (GL) equation: du/dt = u — u 3 + d 2 u/dx 2 [HI]. 
However, it is different from the GL system in that the final state is a uniform state in the GL system and the final 
state is a two-domain state in our system. It is a unique point in our model that the two-domain structure of equal 
size appears naturally for any system size L. 

Equation (3) can be generalized into an asymmetric model: 



dX 1 cXi + d^Xa) 



dt 1 + aXx + bX 2 
dX 2 = cX 2 +d 2 (X 1 ) 
dt 1 + aX 2 + bX x 



-X 1 + D X V 2 X 1 , 

-X 2 + D X S7 2 X 2 , (6) 



where d\ and d 2 are assumed to take different values. Figure 3(a) diplays stationary profiles of X\ (solid curve) and 
X 2 (dashed curve) at a — 0.4, 6 = 5, c = 2, d\ = 0.4, d 2 = 1, and Dx = 1 for L = 100. The sizes of the AVdominant 
region and the ^-dominant region are not the same in this asymmetric model. However, the maximum and the 
minimum values of X\ and X 2 and the width of the domain walls are almost the same, which are expressed as 
X m ax and X min . If d\{X 2 ) and d 2 (Xi) are assumed to takes constant values and d\{X 2 ) < d 2 {Xi), a domain wall 
between a Xi-dominant region and a ^-dominant region moves as the Xi-dominant region shrinks, which is similar 
to the motion of the domain wall in the asymmetric Ginzburg-Landau equation: du/dt = u + eu 2 — u 3 + d 2 u/dx 2 . 
However, if the AVdominant region shrinks, (Xi) decreases and (X 2 ) increases in time in our system. This is due to a 
negative feedback effect involved in our system. Finally, the domain walls become stationary, when d\(X 2 ) — d 2 {X\) 
is satisfied. If the system size L is sufficiently large, the ratio r — l/L of the domain size / of the Xi-dominated region 
is evaluated from the condition d\(X 2 ) = d 2 (Xi) as 

di{(l - r)X max + rX min } = d 2 {rX max + (1 - r)X min }, (7) 

because (X 2 ) ~ (1 — r)X max + rX m i n and (Xi) ~ rX rnax + (1 — r)X m i n , if the width of domain walls is not taken 
into consideration. If X m i n X max , r is approximated at r = d\/{d\ +d 2 ). Thus, the size ratio is determined by the 
ratio of the parameters d\ and d 2 and does not depend on the system size L. We studied the control of the domain 
size in the Ginzburg-Landau type equation in the previous work, which was the same mechanism as the present 
one fl6l ]. Figure 3(b) displays numerically obtained size ratio r's as a function of di at a = 0.4,6 = 5,c = 2,d 2 = 1, 
and Dx = 1 for L = 100. The dashed curve denotes r = l/(e?i + 1). Fairly good agreement is seen. The ratio 
control can be attained by the choice of the parameters d\ and d 2 . This is a mechanism of the ratio control of our 
differentiation model. The feedback effect via the Y- variable with infinitely large diffusivity controls the ratio of the 
different cell types. In this asymmetric model, a two-domain structure appears naturally, and the ratio of the domain 
size is uniquely determined by the system parameters. It might be applicable to the ratio control in the early stage 
of the slime mold. 



III. A CASCADE MODEL OF CELL DIFFERENTIATION 



Next, we construct a new cascade model of the competitive reaction-diffusion equations as shown in Fig. 4(a) based 
on Eqs. (3) and (6), although another cascade model based on a more complicated version of Eq. (1) was proposed 
by Meinhardt to study the hierarchical subdivision into gap-gene, pair-rule gene and segment polarity gene in the 
Drosophila [lCj. In Fig. 4(a), a system of fourteen elements in three layers is shown. Each element denotes a protein 
(morphogen) with number i which determines the ith cell type. The horizontal dashed lines represent the competitive 
interaction between two elements as in the previous model, and the solid lines from the upper layer to the lower layer 
represent an active interaction from a upper- level element to a lower-level element. A negative feedback from the 
lower-level to the upper level might be important, but the feedback effect is not considered in our model for the sake 
of simplicity. The active and competitive interactions are represented by the Hill type equations. For example, a 
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FIG. 4: (a) Hierarchical network of fourteen elements in three layers. The solid lines denote the active interaction and the 
dashed lines denote the competitive interaction, (b) Stationary profiles of Xi (i = 1, • • • , 14) in the cascade model similar to 
Eq. (8) in a one-dimensional space. The system is composed of fourteen elements in three layers. The system size is L — 96. 
The parameter values are a — 0.4, 6 = 5, c = 2, d = 1, a = 0.4 and Dx = 0.1. (c) Stationary profiles of Xi (i = 1, • ■ • , 14) at 
a = 0.4, b = 5, c = 2, di = 0.8, efe = 1, a = 0.4 and D x = 0.1. 



model equation for a system of six elements in two layers is written as 

dX 1 cX^d.jX,) 

~df - l + aX 1 + bX 2 - Xl+DxV Xl > 

dXt = cX a + M Xi) _ X2 + Dx ^x 2 , 
dt l + aX 2 + bX 1 



dX 3 _ aX 1 {cX 3 + d 1 {X A )) 



dt 


l + aX 3 


+ bX 4 


dX± 


all (cA 4 4 


-d 2 (X 3 )) 


dt 


l + aA 4 


+ bX 3 


dX 5 


aX 2 (cX 5 4 


■d!{X 6 )) 


dt 


1 + aX 5 


+ bX 6 


dX 6 


aA 2 (cA 6 4 


-d 2 {X 5 }) 


dt 


1 + aXe 


+ bX 5 



X 3 + D X V 2 X 3 , 



A 4 + DxV 2 X±. 



X 5 + D X V 2 X 5 , 



X 6 + D X V 2 1 6 , (8) 

where a is the coupling constant of the active interaction from the upper layer to the lower layer. Owing to the 
competitive interaction, X 2 is almost zero in the JTi-dominant region. At the domain, X§ and Xq are also almost 
zero because X§ and Xq are activated by X 2 . On the other hand, the Xi-dominant region is separated into the 
JT3-dominant region and the ^-dominant region because of the competitive interaction between X 3 and X±. The 
size ratio of the two regions are also determined by the ratio of d\ and d 2 . We have performed numerical simulations 
of a system of fourteen elements in the three layers. It is a one-dimensional system with system size L = 96. Figure 
4(b) displays stationary profiles of Xi (i = 1, • • • , 14) for the three layer system at a = 0.4, b = 5, c = 2, d = 1, a = 0.4 
and Dx = 0.1. Firstly, the whole space is separated into the two domains: a Xi-dominant region and a X 2 - 
dominant region. Next, the Xi-dominant region is separated into the two domains: a JT3-dominant region and a 
^-dominant region, and the ^-dominant region is separated into the two domains: a As-dominant region and a 
Ag-dominant region. Similarly, the ^-dominant region is separated into Ay-dominant and Ag-dominant regions, the 
A4-dominant region is separated into Xg-dominant and Aio-dominant regions, the As-dominant region is separated 
into An-dominant and Ai2-dominant regions, and, the Ag-dominant region is separated into Ai3-dominant and A14- 
dominant regions As a result, a hierarchical pattern appears. The domain size decreases as 1/2, 1/4 and 1/8 as the 
layer is lowered. This is a kind of binary-tree decomposition of the whole space. If the parameter a is assumed to 
take a different value otk for each layer, and is set to be at = (c + d/2)/{X max (c + d/2 h )} for the fcth layer, peak 
values of A, take almost the same value X max for each layer, because (Aj) = X max /2 k in the A;th layer. For the 
well tuned parameter values of o^, a self-similar binary decomposition of dominant regions will occur. If d\ ^ d 2 , 
a multi-fractal-like pattern appears, because the size ratio is r, 1 — r in the first layer, r 2 ,r(l — r), (1 — r) 2 in the 
second layer, and r 3 , r 2 (l — r), r(l — r) 2 , (1 — r) 3 in the third layer. Figure 4(b) displays a stationary profiles of Xi 
(i = 1, • • • , 14) for d\ = 0.8 and d 2 — 1 at a — 0.4,6 = 5, c = 2, d = 1, a = 0.4 and Dx = 0.1. Inhomogeneous 
decomposition of dominant regions is clearly seen. 

Finally, we have performed a numerical simulation in two dimensions. The system size is 48, and the periodic 
boundary conditions are imposed. The network of the interaction is the same as before in Fig. 4(a), that is, the 
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FIG. 5: Hierarchical differentiation of the fourteen elements with the three layers in a two-dimensional space. The X;-dominant 
regions (i = 1, • • • , 14) are depicted with different colors in a square of size 48 x 48. The parameter values are di = 0.8 and 
c?2 = 1 at a — 0.4,6 = 5, c = 2,di = 0.4,^2 = l,a = 0.4 and Dx = 0.2. A hierarchical structure is clearly seen. The colored 
regions imply that Xi satisfies Xi > 1. (a) Xi-dominant and A2-dominant regions, (b) X3, X4, Xs, and ^-dominant regions, 
(c) Xy, Xs, Xs, Xio, Xn, X12, X13, and Xi4-dominant regions. The number i indicates the cell type. 
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FIG. 6: Time evolution of the number ratio of AVdominant regions satisfying Xi > 1 for i = 1, 2, 3, 4, 7 and 8. 



number of elements are fourteen and the layer number is three. The parameter values are d\ = 0.8 and e?2 = 1 at 
a = 0.4,6 = 5, c = 2,di = 0.4, d,2 = l,a = 0.4 and Dx = 0.2. The initial conditions are random between 0.5 and 
0.6. Figure 5(a), (b) and (c) display patterns at t = 20000 which represent the dominant region satisfying Xi > 1 
for each i using different color for (a) i = 1 and 2, (b) i — 3, • ■ • , 6, and (c) i = 7, • • ■ , 14. These states are final 
stationary states after the transient time evolution as shown later in Fig. 6 and Fig. 7. The whole region is divided 
into two domains: a Xi-dominant region and a ^-dominant region. The Xi-dominant region is a circular blue region 
which is located near the four edges in Fig. 5(a). The area ratio of the Xi-dominant region and the ^-dominant 
region is 0.42 : 0.58, which is close to d\ : d,2 = 0.44 : 0.56. The Xi-dominant region is subdivided into two domains: 
a J^-dominant region and a ^-dominant region. The area ratio of the two regions is 0.44 : 0.56, which is also 
close to the ratio d\ : di — 0.44 : 0.56. Similarly, the X3-dominant region is further subdivided into a ^-dominant 
region and a Xg-dominant region. The area ratio of the two regions is 0.42 : 0.58, which is also close to the ratio 
di : d 2 = 0.44 : 0.56. Thus, the hierarchical differentiation is clearly seen. And it is shown that the area ratio is 
determined by the ratio of d\ and d 2 . 

Figure 6 displays the time evolution of the number ratio of the X^-dominant regions for i = 1,2, 3, 4, 7 and 8. The 
number ratio becomes a stationary values at t ~ 30 for i — 1 and 2, at t ~ 70 for i = 3 and 4, at t ~ 1100 for i = 7 
and 8. That is, the stationary state is attained sequentially from the upper layer to the lower layer. Figures 7(a), 
(b) and (c) show three snapshot patterns of the X%-,Xi- : A5- and ATg-dominant regions at t = 50, 100 and 2000. The 
dominant regions are initially small and random, and the size and the location are determined gradually. Figure 7(c) 
is almost similar to Fig. 5(b), which means that the time evolution is almost stationary at t = 2000. 
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(a) (b) (c) 




FIG. 7: Three snapshot patterns of the Xi dominant region satisfying Xi > 1 for i = 3.4.5 and 6 at (a) t = 50, (b) t — 100, 
and (c) t = 2000. The snapshot at t = 20000 is shown in Fig. 5(b). The number i indicates the cell type. 



IV. SUMMARY 



Competitive interactions and hierarchical interactions are typical in gene networks. We have proposed a simple 
cascade model of competitive reaction-diffusion equations for the cell differentiation. We have first reconfirmed a 
spatially-periodic pattern in a simple one-dimensional competitive reaction diffusion equations, which was originally 
proposed by Meinhardt and Gierer. It is a kind of the Turing pattern, but the origin of the formation of the periodic 
pattern is the long-range enhancement of the competitor. If the diffusion constant for Y- variable is assumed to infinity, 
we get a new model, in which a two-domain structure appears naturally from an random initial condition, and the 
size ratio of the two domains can be controlled by system parameters owing to the negative feedback effect of the 
domain size. Next, we have generalized the model to a hierarchical model. The AVdominant regions appear in a 
cascade manner from the upper layer to the lower layer. The ratio of the AVdominant regions are well controlled by 
the system parameters. Our simple cascade model is not a realistic model based on biological experiments, however, 
it is a useful model to consider the Turing pattern, the ratio control, and the hierarchical differentiation in a unified 
manner. We expect that our model might be applicable to specific cell differentiation processes by changing the 
reaction network and modifying the system parameters suitably. 



[1] G. Nicolis and I. Prigogine, Self-Organization in Nonequilibrium Systems (John Wiley and Sons, New York 1977). 
[2] A. Turing, Phil.Trans. R. Soc.London B 327, 37 (1952). 

[3] H. Meinhardt, Models of biological pattern formation (Academic Press, London, 1982). 
[4] J. D. Murray, Mathematical Biology (Springer- Verlag, New- York, 1989). 
[5] S. Kondo and R. Asai, Nature 376, 765 (1995). 
[6] J. T. Bonner and M. K. Slifkin, Am. J. Biol. 36, 727 (1949). 
[7] R. R. Kay, Seminars in Cell and Devel. Biol. 10, 577 (1999). 
[8] Y. Maeda and M. Maeda, Exp. Cell Res. 84, 88 (1974). 
[9] S. F. Gilbert et al., Developmental Biology (Sinauer, 1997). 
[10] H. Meinhardt, J. Cell. Sci. Suppl. 4, 357 (1986). 

[11] E. Mjolsness, D. H. Sharp, J. Reinitz, J. Theor. Biol. 152, 429 (1991). 

[12] H. Meinhardt and A. Gierer, J. Theor. Biol. 85, 429 (1980). 

[13] A. V. Hill, J. Physiol. 40, R4 (1910). 

[14] M. B. Elowitz and S. Leibler, Nature 403, 335 (2000). 

[15] K. Kawasaki and T. Ohta, Physica A 116, 573 (1982). 

[16] H. Sakaguchi, Phys. Rev. E 64, 047101 (2001). 



